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Abstract 

We derive a set of new unidirectional evolution equations for photonic 
nanowires, i.e. waveguides with sub-wavelength core diameter. Contrary 
to previous approaches, our formulation simultaneously takes into account 
both the vector nature of the electromagnetic field and the full variations 
of the effective modal profiles with wavelength. This leads to the discovery 
of new, previously unexplored nonlinear effects which have the potential 
to affect soliton propagation considerably. We specialize our theoretical 
considerations to the case of perfectly circular silica strands in air, and 
we support our analysis with detailed numerical simulations. 

1 Introduction 

Evolution equations that are able to accurately describe the propagation of 
ultra-short pulses in strongly nonlinear media are extremely valuable tools in 
modern computational nonlinear optics, and especially in fiber optics [T]. On 
one hand, it is always possible to directly solve Maxwell's equations in presence 
of linear dispersion and nonlinear polarization 2 , capturing the full complexity 
of the nonlinear propagation, but in the majority of cases very large computa- 
tional resources are needed for this, which means that one can only effectively 
simulate short propagation distances along an optical fiber. On the other hand, 
it is possible to drastically reduce the complexity of Maxwell's equations by 
using the envelopes of the electromagnetic (EM) fields, while still maintaining 
an accurate description of the important dynamical degrees of freedom. This 
can lead to different approximation schemes based on various different versions 
of the nonlinear Schrodinger equation (NLSE), that have been used quite suc- 
cessfully in the past to describe a variety of linear and nonlinear effects in 
large-core optical fibers. Amongst them, we can mention phenomena such as 
soliton propagation [3] , modulational instabilities and four- wave mixing [3] , and 



1 



third and higher-order dispersive effects pQ . A specific extended version of the 
NLSE (also called the generalized NLSE, or GNLSE for brevity), which also 
includes terms describing the Raman effect [5], self-steepening pQ and the full 
complexity of the group velocity dispersion (GVD) of the waveguide, has been 
used in recent times to describe the important phenomenon of supercontinuum 
generation (SCG), i.e. the generation of a broad coherent spectrum from a short 
input pulse, with spectacular success [6]. This equation has led to a number of 
important advances in the theoretical understanding of SCG in ultra-small core 
fibers such as tapered fibers (TFs) [7] and photonic crystal fibers (PCFs) [5112]. 
In particular the phenomena of emission of solitonic dispersive radiation (also 
called resonant radiation or Chcrcnkov radiation [TO]), and the stabilization of 
the Raman self-frequency shift (RSFS) of solitons by means of such radiation 

1 1 2 j , are recent important advances that have been possible only by means 
of a thorough theoretical understanding of the structure of GNLSE. 

Equations like the GNLSE are based on several assumptions and approx- 
imations. The most important approximation used in all NLSE-based evolu- 
tion equations is the slowly-varying envelope approximation (SVEA) in both 
space and time [13]. This amounts to decoupling the fast oscillating spatial 
and temporal frequencies from the envelope function, which is the only complex 
function used to describe the dynamics. This reduction is only possible if the 
envelope profile contains many oscillations, which is typically the case in the 
visible and mid infrared regime when the pulse duration is larger than a few 
tens of femtoseconds. SVEA allows to reduce the second-order wave equation 
for the electric field into a first-order equation (which only contains the forward- 
propagating modes and neglects the backward-propagating ones) in which the 
evolution coordinate is z, the longitudinal spatial coordinate along the fiber. In 
this way, the computational complexity of numerical simulations is drastically 
reduced. The most crucial assumption on which the GNLSE is based is that 
the z-component of the electric field is very small with respect to its transverse 
components (weak guidance regime). This allows to approximately decouple the 
transverse dynamics from the longitudinal dynamics. The weak guidance regime 
represents a very good assumption for large-core fibers (when the core size is 
much larger than the wavelength of light) and for fibers with a low contrast 
between the refractive indices of the core and the cladding. 

Recently, however, new types of waveguides with a sub-wavelength size of 
the core (broadly referred as photonic nanowires) have become accessible to 
fabrication, which both possess a complex cladding structure, that allows one 
to support solid cores with a sub-wavelength diameter, and have a strong con- 
trast between the refractive indices of the core and the cladding. To this class 
of waveguides belong some specific examples of PCFs, such as extruded PCFs 
[13]; the extremely small interstitial features of Kagome hollow-core PCFs [T3] : 
and TFs, i.e. silica rods with sub-micron diameters surrounded by air or vac- 
uum [7]. As mentioned, the common feature of all these structures is the non- 
negligible nature of the longitudinal component of the electric field (strong guid- 
ance regime) [16] . The importance of such component in the dynamics of light 
has not been clearly recognized in the past, until very recently [16( 117] . In 
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these latter works the authors demonstrate that the calculation of the nonlin- 
ear coefficient of photonic nanowires performed in the scalar theory [T8l [T9l [20] 
underestimates its real magnitude of approximately a factor 2, while the cor- 
rect result can be found by using the more complete vector theory [TO]. Such 
interesting conclusion is also supported by recent experiments performed on 
ultra-small core, highly nonlinear tellurite PCFs [2Tj . 

Inspired by these recent works, in the present paper we wish to further ad- 
vance the theoretical understanding of the evolution equations in sub- wavelength 
core structures. Here we derive a new, logically self-consistent forward-evolution 
equation for photonic nanowires, based on SVEA, that simultaneously takes into 
account (i) the vector nature of the EM field (polarization and z-component in- 
cluded) , (ii) the strong dispersion of the nonlincarity inside the spectral body 
of the pulse (i.e. the variation of the effective modal area as a function of fre- 
quency), and (iii) the full variations of the mode profiles with frequency. In all 
relevant previous works on the subject [321 HH \17\ [22l [23[ [24], all these prop- 
erties were not included in the derivation at the same time. Considering (i), 
(ii) and (iii) altogether and with a minimal amount of assumptions, leads to 
new qualitative results that have not been possible to investigate with previous 
formulations. In particular, we demonstrate through a series of accurate sim- 
ulations that (iii) leads to the emergence of new nonlinear terms that are also 
indirectly responsible for appreciable modifications in the dynamics of RSFS, in 
particular a suppression of the latter. Thus, even though the vector nature of 
the EM field in photonic nanowires may lead to an enhancement of the nonlinear 
coefficient in nanowires Hi . our results show that this is actually counterbal- 
anced by terms and effects that partially suppress such enhancement. However, 
our findings are not in contradiction with previous works on the subject, for 
reasons that we shall explain in the text. 

The structure of the paper is the following. In section[2]we detail the deriva- 
tion of our master equation, Eq. (fT9f . In section [3] we describe the symmetries 
of the master equation, and we demonstrate that our model coincides with pre- 
viously published models under various limiting cases. In section 0] we approx- 
imate Eq. (|19[) for long pulses, and find that new nonlinear terms arise, which 
were not considered in any previously published formulation. Such terms be- 
come increasingly important when decreasing the size of the core, and for a given 
core diameter their impact is maximized around a certain optimal wavelength. 
Finally, in section [5] wc show results of numerical simulations, which underline 
the difference between our new theory and other known models. Conclusions 
are given in section [SJ 

2 Derivation of master equation 

In this section we present a detailed derivation of our master equation, Eq. (|19j) . 
Our starting point are Maxwell's equations written for a non-magnetic medium 
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in the Lorentz-Heaviside units system [25] : 

V x H = -d t T> + -dtPNL (1) 

c c 

VxE = --ftH (2) 

c 

V-H = 0, V-D = 0, (3) 

where E(r, t) and H(r, t) are the electric and magnetic field respectively, t 
is time, r is the coordinate vector, c is the speed of light in vacuum and 
Pnl(i*, t) is the nonlinear polarization field that we shall specify later. The 
displacement field is given by D(r,i) = e ® E, where e(r±,t) is the dielec- 
tric function that models the linear response of the medium, i.e. its disper- 
sion. Its Fourier transform e(r±,uj) gives the square of the refractive index 
along the transverse direction rj_. The symbol ® indicates a convolution: 
a ® b = a(t - t')b{t')dt' = b(t - t')a(t')dt' = b <g> a. In the frequency 
domain, one has D(r, w) = e(r±, w)E(r, uS). 

In absence of nonlinear polarization (Pnl = 0), the fiber has linear eigen- 
modes in the form of plane waves, solutions of the linearized Maxwell equations 
at a specified frequency u) [26] . These are given by: 

E^^r.^^le^f^-' (4) 

V ^mui 

and 

H„_(r,<) = ^'e^M^', (5) 

V "mu 

where m is a collective index that indicates both the polarization state and the 
mode number, and /3 m (to) is the propagation constant. N mu) is a mode- and 
frequency-dependent normalization constant that will be defined shortly. Func- 
tions e m>a ,(rx) and h. m ^{r±_) give the profile of the fiber eigenmodes [55]. Note 
that e, h have all three components, including a z-component which becomes 
comparable with the other two transverse components for a sufficiently small 
core and for a relatively large contrast between the core and the cladding re- 
fractive indices, as pointed out in Ref. [16]. Fig. [T] shows contour plots of the 
transverse and the longitudinal components of the electric field for two circu- 
lar silica strands of diameters d = 2 /jm [Fig. (Ha,b)] and d = 0.6 /jm [Fig. 
[]Jc,d)] respectively, for a pump wavelength A = 1 /im. Note the gradual ex- 
pulsion of the total (transverse and longitudinal) field from the core towards 
the low-index medium (air in this case). For PCFs and for fibers with a highly 
non-trivial cladding structure, the vector profiles of the electric and magnetic 
eigenmodes must be calculated numerically by means of a mode solver, but for 
the specific case of optical fibers with a perfectly circular core, e and h are 
known analytically, provided that the propagation constant is known, see Ref. 

The full fields are expanded in terms of the eigenmodes as follows: 

E M) = r £ / d^A^^de^^ (6) 
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Figure 1: (Color online) (a,b) Respectively transverse and longitudinal components 
of electric field in a silica strand in air of diameter d = 2 /im, for a pump wavelength 
A = 1 /im. (c,d) The same as (a,b), but for diameter d = 0.6 /jm. The silica core is 
indicated with a dashed black line. Note the expulsion of the total field from the core 
towards the low-index medium. Also, the magnitude of the longitudinal component 
in case (c,d) is more than 30% of the transverse component. 

H(r,i) = f £ / ^W^^e^t")^, ( 7 ) 

where we have used the convention f(t) = [27r]^ 1 J dw/(w)e _i "' for the direct 
Fourier transform and f{ui) = J dtf(t)e lult for its inverse. It is customary to 
define the normalized field e muJ = e mu) / V N mul (and similarly for the magnetic 
field), in order to have an orthonormal basis on which to expand any physical 
field, analogously to what is done in quantum mechanics. We further note that 
in Ref. [16] the fields are expanded by assuming that the modes e muJ and h mw 
are unchanged around a reference frequency loq, which is arbitrary and taken 
to be equal to the central pulse frequency. Thus in Ref. [16] it is assumed that 
all spectral components inside the pulse will have the same mode profile, and 
there is no integration over the frequencies analogous to Eqs. (UK]), which is 
instead considered in Refs. [23} HH [22]. In the present work we wish to be able 
to describe the full variations of the mode profiles with frequency because, as it 
was already demonstrated in Ref. [27], this actually leads to a large dispersion 
of the nonlinearity that strongly affects the dynamics of the resulting equations, 
and, most importantly, new nonlinear effects that we shall analyze in detain in 
section [4] 

The correct orthogonality relations between the linear eigenmodes are given 
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Figure 2: (Color online) (a) Plot of the dimensionless coefficient r = 71/70 = 
pioocypoooo ag a f unc ti on Q f wavelength, for silica strands in air of diameters d = 
0.5,0.6,0.8,1,1.6 and 4 /j,m in the fundamental mode of one polarization (m = 1). 
This function is predominantly positive in a broad range of wavelengths, r quantifies 
the relative importance between the new nonlinear term T3 in Eq. (|28p and the Kerr 
effect term Tl. (b) Plot of the dimensionless coefficient a s = 4r 1000 /[r 0000 ^oTR] = 
471 / ' [jo^oTr] as a function of wavelength, as quantifies the relative importance be- 
tween the new nonlinear term T3 in Eq. (|28p and the Raman effect term T2. Both r 
and as tend to zero at all wavelengths when progressively increasing the core diameter. 



by the following cross product integral [26) : 

i J [e ju (r ± ) x h^(r ± )] ■ zdr ± = i / [e* ku (r ± ) x h ju (r x )] ■ zdr x = 5 jk N juJ 

(8) 

where Nj U is a normalization constant. However, what is ultimately important 
is not the absolute normalization, but the relative change of the normalization 
constant (and thus of the modal effective area) with wavelength. That is why 
we decide to use in the following, without loss of generality, the orthonormal 
quantities e, h, which satisfy the simpler conditions 



e ju (r±) x h^(r x ) 



• zdr 1 = S 



(9) 



It is customary jTSJ [23] to define a generalized Poynting vector as follows: 
S mu {r,t) = -e c[E mw (r,t) x H*(r,t) +E*(r,i) xH mw (r,t)]. (10) 



By using Maxwell's equations (p~|[3]) and the vector theorem V • (a x b) = b 
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(V x a) — a • (V x b), one can find the divergence of S: 

V-S mw = -eo {iwH* ■ H mw — E maj • <9 t D* — E ma) • dtP^L — H TOW • 9 t H* + iuie^E* ■ E TOW } 

(11) 

After integrating both members of Eq. pip in the transverse dimensions, per- 
forming a time integration as in (23] , and by using the reciprocity theorem |26j 
one obtains the following forward-propagating equation (first derived by Kolesik 
et al. in 2004, see Refs. [23]): 

d z A mul = ^V^- 2 J e* mu ■ P NL , w (r)dr ± . (12) 

Equation (j 1 2|) is the basic equation for the common Fourier component of the 
EM fields, on which any forward-propagating theory must be based. Note that 
Eq. (|12[) is written for a generic frequency uj, thus no arbitrary reference fre- 
quency has been selected up to this point. Attempts to directly solve Eq. (fT2"j) 
have been given, in various forms, in Refs. [33] and [53]. However, this approach, 
although exact in absence of backward-propagating modes in the waveguide, 
does not provide a transparent understanding of the physical processes involved, 
and it does not lead to any serious reduction of the computational costs with 
respect to solving directly Maxwell's equations. 

To further progress from Eq. (|12|) . one needs to write explicitly the nonlinear 
polarization field for silica in terms of the electric field components. This can 
be done in several ways, depending on which particular medium or model one 
is interested in. For instance, in Ref. [TB] the authors consider the vector Kerr 
nonlincarity of silica in absence of Raman effect. In Ref. [T7] the Raman effect 
has been included, but emphasis was given to the analysis of stimulated Raman 
scattering. For our purposes we choose the following form for the nonlinear 
polarization of silica (see also Ref. [24]): 

Pnl = \x£L x (r±){v(t) J R(t-t 1 )[E(t 1 )-E*(t 1 )]dt 1 + J R(t-t 1 )E(t 1 ){V(t)-V*(t 1 )]dt 1 + 

+ f R(t-t 1 )E*(t 1 )[E(t)--E(ti))dh}, (13) 

where Xxxxx (r_i_ ) is the third-order susceptibility coefficient, which is the only 
approximately independent component of the susceptibility tensor that survives 
the symmetries |28j . and is a function of the transverse coordinates. R{t) is the 
Raman-Kerr convolution function that describes both the instantaneous (Kerr) 
and non-instantaneous (Raman) nonlinearities in the medium. In the case of 
silica one has: 

R(t) = (l-6)5(t) + 9Q(t)h(t), (14) 

where h(t) = [rf + t|][ t i t 2] 1 ex P(~V r 2) sin(i/ri) is the Raman response func- 
tion set by the vibrations of silica molecules induced by the EM field, with 
T\ = 12.2 fs and — 32 fs. Coefficient 9 parameterizes the relative importance 
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between Raman and Kerr effect in Eq. (| 14[) . and the experimental value for 
silica is 9 ~ 0.18 pQ. In Eq. (TT4|) , <5(i) is the Dirac function and 9(i) is a 
Heaviside (step) function that is introduced to preserve causality. In addition, 
the response function is normalized in such a way that f_™ R(t)dt = 1. Other, 
more precise forms of the Raman response functions for silica are known pQ. 
In a more accurate formulation, the isotropic and the anisotropic parts of the 
nuclear response should be included [29]. However, the following discussions are 
largely independent of the precise form of R(t), and are valid in general for any 
functional form of the response function. 

After substituting expression (TT5|) into Eq. (TT2"|) , we explicitly obtain: 



OA -— - 1 



/ X ( *Lxdr±(j2 I dJdJ'R{u-J)A k ^A l ^„A 

kdl 



Ac 4 (2tt) 2 

' • (15) 



Only at this point we choose an arbitrary reference frequency lu ; furthermore, 
in order to capture the variations in the mode profile around ujq, we expand the 
linear modes e in Taylor series: 



Alu 

LOO 



(16) 

where Alu = cu — lu is the frequency detuning from the reference frequency luq, 

( 7*1 

and the quantity fm,w Q is given by the term inside the square brackets in the 
last member of Eq. (fTB]). and is proportional to the j-th derivative of the mode 
profile. A Taylor expansion equivalent to our Eq. (|16p was also proposed in 
Eqs. (28) and (37) of Ref. [53]. However, in this latter work the author (being 
mostly interested in large core, low-contrast PCFs) gets rid of the longitudinal 
components of the EM fields [see his Eq. (13)], analyzes the scalar case, and 
does not elaborate further on the physical meaning of the corrections induced by 
the higher-order terms. Expansion (fTB"]) is also not considered in Refs. [TBI [30] , 
In our derivation, the expansion (fT6|) is a crucial step. In fact, by taking into 
account the modifications of the quasi-linear mode profiles, we are led to a 
correct identification of new nonlinear effects that, under certain circumstances, 
strongly compete with the RSFS. We shall describe the latter effect in section 
® 

By inserting expression (116|) into Eq. (|15|) . and after integrating over the 
transverse coordinates of the fiber, we can single out the following constant 
nonlinear coefficients: 

frf*(j) fO)irf(p) f*(«)i I [f*(i) f(p)irfO) f*(«)i i rf*(i) e*( v )]\c( h ) f(p)]\ 

irm 
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We shall analyze this complicated but very important object and its large 
symmetries in the next section. T^^ d [in the following, the upper indices jhpv 
shall always indicate the order of the derivatives in the Taylor expansion in 
(fTB|) , while the lower indices mkld shall always denote the polarization and 
mode states] contains information about the nonlinearity experienced by the 
different modes or polarizations in the fiber. Note that the structure of T^j^ 
closely resembles the one of the nonlinear polarization Pnl inside Eq. (fTS)) . 
In addition, all components of r^ 1 ^ have the same physical dimensions. In 
particular, it is convenient to use the usual units for nonlinear coefficients of 
W^mT 1 for the quantities k T J ^i d . 

The final step in our derivation is to find an equation in the time domain 
for the envelope of the electric field. In order to do this, we decouple the fast 
oscillating phase by introducing the envelope function Q, therefore using SVEA: 

A mu = Q ro , Aw e^o-/^> e -^ot. (18) 

In this paper we shall deal with perfectly circular photonic nanowires, i.e. silica 
strands surrounded by air or vacuum. In order to simplify our treatment we 
assume that the waveguide is operated in a single-mode regime, so that the index 
m = 1, 2 only represents the polarization states of the fundamental mode, which 
in perfectly circular fibers are degenerate \fi\{uS) = /^(w)] and orthogonal in the 
sense of the cross product, Eq. (|9|). This assumption can of course be relaxed 
without problems in our formulation, but here we wish to focus on the novel 
properties of our equation. After some regrouping and a Fourier transformation 
over the new independent variable Auj, wc finally obtain our master equation: 

id z Q m + D m {id t )Q m + fco E E r S [G J W ® *^(t)] = 0, (19) 

kid jhpv 

where fco = ljq/c is the vacuum wavenumber, D m (idt) = /3 m (^o + idt) — Pm{<^o) 
is the dispersion operator that encodes all the complexity of the fiber GVD 
around the reference frequency [23j . and 

is a function that naturally contains the dynamics of the shock term [through 
the term (1 + Aui/ojq)]. The convoluted nonlinear fields in Eq. (JTUJ) are given 
by: 

^hwM _ [(idt) h Qk(t)] {R(t) g {[{idt) p Qi{t)} \{-id t YQ* d {t)])} , 01 v 

^kld W — h+p+v ' \ lL > 
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that also depends on polarization state indices and derivative orders in the 
expansion of Eq. p^|) . 

Equation (fl9]) is the central result of the present work. It differs from other 
formulations in that it has been derived by making no assumptions or approx- 
imations on the frequency dependence of the field profiles of the linear modes 
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of the waveguide. Higher-order derivatives of the mode profiles embedded into 
Eq. (|19p are responsible for the appearance of additional nonlinear terms in the 
equations. These new terms obviously become increasingly small with increas- 
ing orders of the derivatives, since for every derivation there is a factor 1/ljq. 
However, as we demonstrate in section^ the new first-order term has a strength 
comparable with or sometimes even larger than the Raman and Kerr terms, and 
there are circumstances in which it must have observable consequences. 



3 Symmetries of the master equation and reduc- 
tion to known models 

The nonlinear coefficient T^ 1 ^ given in Eq. (fTF)) possesses very large symme- 
tries that we wish to analyze in this section, following similar arguments as 
in Ref. |30j . Some of these symmetries are intrinsic, i.e. valid for any fiber 
geometry, while some other are strongly geometry-dependent. As we shall see, 
the symmetries of T^ 1 ^ can be effectively used to reduce enormously the com- 
plexity of Eq. (fl9)) . As a practical important example, here we specialize the 
description of such symmetry relations for the case of a silica nanowire with 
a perfectly circular core embedded in air. In this case, we consider only the 
fundamental mode of the fiber, which has two polarizations m = 1,2 with the 
same propagation constant = fci^)]. In the following we always assume 

that when index m is fixed, then n = 3 — m denotes the opposite polarization. 
As before, indices j, h,p, v > are used to denote the orders of the derivatives 
in the definition (TTT)) of T, and they can take any positive integer value. 

Directly from the definition (jTTJ) , one immediately gets the relation r^ 1 ^ = 
r*?;Jj) so the two lower an upper central indices can be simultaneously inter- 
changed without changing the value of the nonlinear coefficient. Other sim- 
ple relations (that can be directly verified by the reader) are the following: 

pjhpv _ yvphj jyjhhh _ yhjhh yjhhh _ yhjhh All t hp abovp enuali- 

mklm mlkm^ mmnn mmnni mmmm mmmm' o4uc*n 

ties are always valid for any geometry of the waveguide. 

In the special case of perfectly circular fiber, from direct computation of 

Eq. (p~7]) by using the explicit formulas for the two (even and odd) degenerate 

polarizations of the fundamental mode given in Ref. , one obtains the useful 

enualities V jhpv - V jhpv V jhpv - V vhpj and rihhh _ yhjhh 

equalises 1 mkld ~ 1 (3-m)(3-fc)(3-i)(3-d) ' 1 mkld ~ 1 dklm' duu 1 mnnm ~ 1 mnnm- 

Most importantly, many of these coefficients vanish identically: T^ p ^ nm = 

T^jhpv j^jhpv yjhpv j^jhpv T^jhpv ~njhpv T^jhpv n 

mmnm mmmn nmmm nmnn nnmn nnn m mnnn 

We now demonstrate that our propagation equation (TIT)]) can be reduced to 
the previously published model of Ref. [16] by using the above symmetries, and 
by taking into account only the lowest order in the modal profile expansion Eq. 
(|16[) . For this purpose, let us write the zero-th order approximation (j = h = 
p = v = Q) of Eq. (H9|): 

id z Q m + D(id t )Q m + fc £ r« [G° $ooo] = o. (22 ) 

kid 
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We now explicit the summations in Eq. for the nonlinear functions 
obtaining: 

id z Q m + D(id t )Q m + k r™™ mm [g° ® {Q m (jz® |Q m | 2 )}] + fc r™ ro [g° ® {Q„ (#® |Q m | 2 )}] + 

fcor™^nm [G° ® {Qm (R ® Qn<9*„)}] + fcoCZm [G° ® {Q m (i? ® QmQ£)}] + 

fcoC°L [g° ® {g„ (n® |q„| 2 )}] + fcoC°°n„ [g° ® {Q™ (fl® IQ„I 2 )}] + 

fc r° m Zn [G° ® {Qn (R <8> Qm<&)}] + fcolTZm [G° ® {Q„ (i2 ® QnQm)}] £») 

According to the symmetries reported above, the fourth, fifth, sixth and seventh 
term of Eq. (|23|) vanish identically. Moreover, the eighth and ninth term have 
an identical nonlinear coefficient. A careful comparison between all terms of Eq. 
([23)) shows that, indeed, the structure of this equation is identical to the one 
reported in Eq. (45) of Ref. [16 , in the limiting case of a single mode regime, 
and without the Raman effect. However, our definition of the envelope Q as 
given by Eq. (|18[) automatically incorporates into the equations the frequency 
variations of the normalization constant (which is proportional to the effective 
mode area, see Ref. [27]), something that has not been done in the approach of 
Ref. [16], where the constancy of N mu around loq was assumed. 



4 New nonlinear effects in photonic nanowires 



The higher-order derivatives in the Taylor expansion performed in Eq. (|16|) 
are expected to produce additional nonlinear terms, previously unknown, that 
become increasingly important when the core diameter is reduced and/or the 
refractive index contrast between cladding and core is increased. Here we at- 
tempt to obtain a physical understanding of these new terms, at least at the 
first-order in the Taylor expansion. In the following we just want to isolate 
and understand as clearly as possible the new nonlinear terms that come out of 
Eq. (flU)) . In order to do this, we just take one polarization state (m = 1) for 
simplicity, we expand Eq. ([TO|) and collect all the zero-th and first-order terms 
coming from the products between the G, T and $ functions, obtaining 



i 

LJQ 



J_a$ooo , $ 100 , $ 010 + $ 001 



id z Q+D{id t )Q+kvT ww & ,m + — k T^d t ^ v +k T 

UJo 

(24) 

By using the relation (i/w )<9 t $ 000 = $ 100 + $ 010 - $ 001 (that can be derived 
just by explicitly calculating the derivative of $ 000 ), one can slightly simplify 
Eq. into 



id z Q + D(id t )Q + 7o 



1 + - 1 + 2^)9, 

\ 7o 



$ uuu + 2 7l $ uui ~0, (25) 



where we have defined the constants 70 = fecT 0000 and 71 = fcoT 1000 . A plot 
of the ratio r = r 1000 /! 10000 — 71/70 is given in Fig. (DJa) for silica strands 
of various diameters. Interestingly, for a given diameter, r is predominantly 
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Figure 3: (Color online) (a) Plots of calculated values of a?i = 47i/(woio) in [m _1 W _1 ] 
describing the dispersion of the first-order nonlinear coefficient as a function of wave- 
length, for silica strands in air of diameters d — 0.5,0.6,0.8, 1, 1.6 and 4 /im in the 
fundamental mode of one polarization (m = 1). Pulse duration is to — 100 fs. For a 
given diameter d, the largest value of |ai| is located at approximately A ~ Cd, with 
C ~ 1.6. (b) Plots of calculated GVD for the same diameters as in (a). Red dashed 
line indicates the GVD for bulk silica. A comparison between (a) and (b) shows that, 
for silica strands, the maximum of |qi| is always located in the normal GVD region. 



positive, and becomes negative only in regions of very short wavelengths, not 
easily accessible experimentally (typically below 200 nm). Also, r is significantly 
larger than unity for wavelengths A > Cd, where C ~ 1.6 for perfectly circular 
silica strands in air. In order to gain some insight into the physical meaning 
of all nonlinear terms in Eq. ([23)) . we consider relatively long pulses (of width 
t ^> 100 fs) for which the envelope evolves slowly along the fiber. Under this 
assumption, one can approximate the convolutions with nonlinear combinations 
of the field Q and its derivatives PQ. In particular, one has 

/+oo r+co 
R{t')\Q(z,t- t')\ 2 dt' ~ Q{z,t) / R{t') {\Q{z,t)\ 2 - t'd t \Q(z,t)\ 2 } dt' = 
- oo J — CO 

= \Q\ 2 Q - T R Qd t \Q\ 2 , (26) 
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+00 



— Q(z,t) / R(t - t')Q(z,t')d t >Q*(z,t')dt' 

LUq 



Q(«,i) / R{t')Q{z,t-t') 



dQ*{z,r) 
dr 



dt' 



Q(z, t) / R(t') [Q(z, t) - t'd t Q(t)] 

LUq 



dQ*(t) f! d 2 Q*(t) 



dt 



dt 2 



dt' 



~ [Q 2 d t Q* - T R Q 2 8 2 Q* - T R Q\d t Q\ 2 + T R Q(d t Q){d 2 Q*)] , 

where T R = f_°° t'R(t')dt' is the first moment of the Raman nonlinear response 

function, and T' R = f_°° t' 2 R(t')dt' is the second moment. The above simple 
Lorentzian model for the Raman response function gives T R ~ 1.5 fs and T' R ~ 
—23 fs 2 for silica fibers. Substituting Eqs. (|26H27p into Eq. (j2"5")l , and neglecting 
the second order derivatives and all terms of order of loq 2 or Tr/loq, after some 
algebra we obtain our final expression (in the first-order of the perturbation 
theory) for the approximate evolution equation for long pulses: 



id z Q + D(id t )Q + 70 




T R Qd t \Qf 



+ ^\Q\ 2 d t Q 



T3 



(28) 

In Eq. (|28|) . the term indicated by TO represents the well-known shock operator 
P] . Terms Tl and T2 in Eq. (|28|) correspond to the usual and easily recognizable 
Kerr effect and RSFS respectively [U [31] . The additional term T3 is connected 
to a nonlinear change of the group velocity, proportional to the field intensity 
\Q\ 2 ■ Such nonlinear term arises quite unexpectedly, and it is strictly related 
to the way in which the mode profile changes with frequency. Slightly different 
frequencies have different linear mode profiles, which in turn will have slightly 
different group velocities. This change becomes nonlinear, because the time 
derivative inside <fi 001 in Eq. (|2"Tj) is embedded into the full nonlinear convolution 
of Eq. (f2"Tj) . thus giving rise to term T3 in Eq. (f2"8"|) . It is now easy to see 
that ratio r of Fig. [2ja) quantifies the relative importance between the new 
nonlincarity T3 and the Kerr nonlinearity Tl in Eq. (|28]) . Note that r can be as 
large as 6 in silica strands inside the frequency window accessible experimentally. 

The relative importance between the T3 and T2 terms is regulated by a 
dimensionless parameter as = 47i/(7oWoTr), which is shown in Fig. [H^b) as 
a function of wavelength for silica strands of various diameters. Again, this 
function is positive above 200 nm. Note that as can be as large as 15 inside 
the frequency window accessible experimentally. 

The presence of the extra first-order term T3 has not been previously re- 
ported in the context of unidirectional pulse propagation equations, and it con- 
stitutes an important result of the present work. It can be demonstrated [32j 
that this term may lead to either a suppression or an enhancement of the RSFS, 
depending on whether the coefficient as is respectively positive or negative. In 



13 



addition, the discussions above demonstrate that there are broad regions of 
wavelengths in which T3 can be several times larger than both the Kerr and the 
Raman terms. Fig. [3]Ja) shows the calculated value of the new first-order non- 
linear coefficient a\ = ^)i/{u>oto) as a function of wavelength, for silica strands 
in air of various core sizes and for to = 100 fs. It is interesting to see that the op- 
timal wavelength for which |ai| is largest approximately corresponds to A ~ Cd, 
again with C ~ 1.6. Moreover, it can be seen in Fig. [3ja) that such a linear 
law seems to be a universal behavior, once that geometry and core and cladding 
materials are specified. Fig. \3[b) shows the calculated GVD of silica strands 
in air of various core sizes. By comparing Fig. OJb) with Fig. [3fa) one can 
see that the condition A ~ Cd is always satisfied in correspondence of normal 
GVD for this kind of circular geometry. In our extensive numerical simulations 
we have seen that the new first-order nonlinear terms make a substantial dif- 
ference in the pulse propagation only when the input pulse is launched in the 
anomalous dispersion, i.e. in presence of soliton formation, where it can en- 
ter in competition with the RSFS of solitons, while in the normal dispersion 
regime the contribution of the term proportional to 001 on the dynamics is 
largely negligible. For this reason, it seems that silica circular nanowires are 
not optimal for evidencing the effects of the new nonlinear term T3 regulated 
by nonlinear coefficient a,\. Possibly this is why such new effects have never pre- 
viously been noticed experimentally, although techniques to fabricate photonic 
nanowires have already been available since a few years [TH] • A detailed study of 
the optimization of the new nonlinear terms in sub-wavelength-core PCFs and 
high-refractive index glasses will be carried out in a separate publication |32j . 

As a final note, we would like to stress again that the expansion of Eq. (I27| 
is valid only for relatively long pulses. For ultra-short pulses, such as to < 50 
fs, all the physical effects and terms that are clearly visible in the formulation 
of Eq. (|27|) will be naturally embedded into the full convolution contained in 
$ 001 , see Eq. flM). 

5 Numerical analysis and comparison with pre- 
vious models 

In this section we show results of the extensive numerical simulations that was 
performed by using Eq. (fT^J). We shall compare the output spectra obtained by 
pulse propagation in different dispersive regimes, and for the situation in which 
Eq. (frnj) is truncated at the zero-th order, as opposed as when instead all terms 
of first-order are included. 

As our representative case, we take a silica strand of core diameter d = 0.6 
[im in air, the dispersion of which is displayed in Fig. \S[.b). For this diameter, 
the two zero GVD points are approximately located at Ai = 0.468 /im and A2 = 
0.718 /im. The input pulse duration is always taken to have a temporal width 
equal to to = 100 fs, with the shape of a hyperbolic secant, and with an input 
power P such that the 'number of solitons' parameter N = [joPto/lPzi^o)]} 1 ^ 2 
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is always equal to N — 7, where /3a (Ao) is the GVD coefficient calculated at the 
pump wavelength Ao- Propagation distances arc given in dimcnsionlcss units: 
£ = z/Ld 2 (Ao), where Lr>2(Xo) = ^/l/MAo)! is the second-order dispersion 
length at the pump wavelength Ao HJ. Also, for simplicity we assume for the 
moment that only one polarization state is excited (m = 1), and we neglect the 
dynamics of the other polarization state. In the final part of this section we shall 
also show some result obtained when taking into account the full polarization 
structure of Eq. (fT!)]) . 
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Figure 4: (Color online) Simulation of light propagation in a photonic nanowire of 
diameter d — 0.6 /im, for only one polarization state (m = 1). Input pulse duration 
is to = 100 fs, and input pulse power is always N = 7. (a) Output spectra of the 
simulated Eq. (|19p in normal dispersion regime after a propagation length of £ = 1, 
by using different truncations of the Taylor series contained into Eq. ((28} . Input pulse 
is centered at Ao = 1.2 /im (black solid line). Blue solid line indicates the solution 
of Eq. I|19p when truncating the sum up to the first-order in the Taylor expansion. 
Green dashed line gives the same as the blue solid line, but when the sum in Eq. (I19|l 
is truncated at the zero-th order of expansion, i.e. when neglecting the influence of 
the new first-order nonlinear terms described in section f5] (b) Same as in (a), but 
for a pump wavelength Ao — 0.8 /mi. (c) Same as (a), but pumping at a wavelength 
Ao = 0.55 /im in the anomalous dispersion, for £ = 1. (d) Same as (c), but for £ = 1.4. 



Fig. SJa) shows the output spectrum of an input pulse (displayed in black 
solid line) propagating in the normal dispersion regime of the nanowire, when 
pumping at Ao = 1.2 /im, and after a dimcnsionlcss propagation distance of 
£ = 1 (Ld2 — 6 mm). At this wavelength we have r = 1.351 and as — 2.295, 
see Fig. r2Ja,b). The blue solid line in Fig. 0|a) shows the result of the simulation 
of equation (|19p , when considering the expansion of the Taylor series of the mode 
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profiles up to the first order. The green dashed line in Fig. 0Ja) shows the results 
of the simulation of Eq. (|19[) . when considering only the zero-th order in the 
expansion. It is possible to observe that in the case of normal dispersion there 
is indeed only a small difference between the green dashed line and the blue 
solid line, which means that the new first-order nonlinear terms contained in 
Eq. (|19p give only a relatively small contribution in this regime, as anticipated 
before. 

When launching the input pulse at Ao = 0.8 /im {Ld2 — 10 cm), still in the 
normal dispersion regime but closer to one of the zero GVD points, all three 
curves nearly coincide, see Fig. EJb). In fact, the new nonlinear coefficient is 
smaller than in the previous case, r = 0.3973 and as — 0.45. Thus we can 
conclude that in absence of soliton dynamics, the two models give qualitatively 
and quantitatively very similar results. 

A completely different scenario is observed in presence of solitons, see Fig. 
|4jc) . This figure shows the output spectrum of a pulse launched in the anoma- 
lous dispersion regime at Ao = 0.55 /im, for a propagation of £ = 1 (L_D2 — 30 
cm). At this wavelength, r — 0.1537 and as = 0.1196. Fig. 01c) shows that 
in the initial stages of the propagation, the stabilization of RSFS by means of 
resonant radiation (at around 800 nm in the figure) emitted near the zero GVD 
point |llj . occurs earlier when the equations do not take into account the new 
first-order nonlinear term T3 in Eq. I|28p . as the strongest soliton red-shifts at 
a faster rate in this latter (less precise) model. It was mentioned in section 0] 
that the novel nonlinear terms associated to the convolution </> 001 in Eq. (| 19[) 
can give rise either to a suppression or an enhancement of the RSFS of solitons, 
depending on whether the sign of coefficient as is positive or negative, respec- 
tively. This is due to a peculiar interplay between the Raman effect and the 
nonlinear change of group velocity induced by the new terms in the solitonic 
regime [32 . Despite the smallness of coefficients r and as at this wavelength, 
a quite strong suppression scenario is already visible in Fig. Hfc). In fact, one 
can see that in the full model that takes into account the additional first-order 
nonlinear terms (blue solid line) solitons arc pushed towards longer wavelengths 
at a much slower rate than in the model without first-order nonlinear effects 
(green dashed line). If the pump wavelength is located in a frequency region 
where as is negative, exactly the opposite scenario would be observed, and the 
blue solid line would be shifted towards longer wavelengths with respect to the 
green dashed line, which means that the RSFS rate is increased. However, such 
a region of positive as is located in the UV in circular silica strands [see Fig. 
Ht[b)] and it is not reachable experimentally. Thus, for the geometry consid- 
ered in this paper, the suppression of RSFS is the only physically meaningful 
scenario. Due to the relative smallness of \oig\ inside the anomalous dispersion 
regime of the circular silica waveguides under consideration, the corrections to 
RSFS induced by the new nonlinear terms are also not as large as they could 
be. However, as we shall show in detail in a future publication [32] . by using 
properly designed PCFs and high-refractive index materials, one is able to bring 
the wavelength for which |ai| has a maximum [see Fig. [3ja)] inside the anoma- 
lous dispersion regime. In such a case, the new nonlinear terms will completely 
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change the dynamics of solitons, leading to qualitatively new regimes of nonlin- 
ear propagation in which the RSFS of solitons can be completely stopped for 
many dispersion lengths, without the use of stabilization techniques based on 
solitonic resonant dispersive wave emission [TT] . 

In the same conditions as in Fig. BJc), but for a slightly longer propagation 
distance £ = 1.4 [see Fig. Bid)], the solid blue and the dashed green lines become 
qualitatively very similar, apart from the visible differences in the dynamics of 
generation of resonant radiation, around 800 nm in Fig. H^d), the bandwidth of 
which is somewhat reduced in the model truncated at the first-order (blue solid 
line) than in the model truncated at the zeroth-order (green dashed line). 

In all the above calculations we have simulated only one of the two polariza- 
tion states, in order to clearly display the differences between several models and 
approximations. However, Eq. (|19[) has a vector nature and actually represents 
a pair of coupled equations for the two polarization fields Q\ and Q2- In Fig. 
[5ja) the final spectrum of a linearly polarized input pulse (95% on the m = 1 
axis and 5% along the m = 2 axis) is shown, simulated by using the full Eq. 
(fl^) . All other parameters for the pulse and the fiber are the same as in Fig. 
|U(c) . Fig. OJb) shows the time domain picture which corresponds to Fig [5ja) . 
As expected, the formation of several vector solitons of different amplitudes is 
observed, each of them subject to RSFS according to their individual intensities. 

6 Conclusions and future work 

In conclusion, we have formulated a new set of vector unidirectional evolution 
equations [Eq. p^|) ] for silica photonic nanowires, i.e. for waveguides with a 
core smaller that the wavelength of light, and with a large step in the refractive 
indices of the core and the cladding, for which the z-component of the electric 
field is not negligible. Our model differs from previous formulations in that 
we simultaneously take into account the vector nature of the EM field and the 
variations of the linear mode profiles with frequency. Such variations generate 
new higher-order nonlinear terms, given by the Taylor expansion in Eq. (|19|) . An 
analysis of the long-pulse limit of Eq. (fTO)) has been carried out, which makes 
the physical understanding of the new first-order nonlinear terms implicitly 
contained in Eq. (fl9| much more transparent. This analysis led to Eq. ([28]) . 
in which a nonlinear group velocity term appears, the strength of which is 
regulated by coefficient a%, which for silica strands in air of given diameter 
d takes appreciable values in the normal dispersion for wavelengths A > Cd, 
with C ~ 1.6. Numerical simulations of pulse propagation in silica strands 
were provided. The difference between the output spectra of our model and 
previous formulations is significant, especially when the propagation occurs in 
the anomalous dispersion regime, where soliton dynamics is mostly affected by 
the new nonlinearities. 

Future works include the use of ultrasmall solid-core PCFs and interstitial 
features in the cladding of hollow-core PCFs, as well as high-refractive index 
nanowires, for the optimization of the novel nonlinear terms. This optimization, 
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which makes use of the great flexibility in the dispersion tailoring of the above 
fibers, will bring the wavelength range at which the new nonlincarities are most 
effective (i.e. where \ a\ | is around its maximum) inside the anomalous dispersion 
regime of the fiber. This will create a truly qualitatively new and unexplored 
nonlinear behavior in such waveguides. 

This work is supported by the German Max Planck Society for the Advance- 
ment of Science (MPG). The authors would like to acknowledge Prof. P. St. J. 
Russell for stimulating discussions, and S. Stark for helping with the numerical 
codes. 
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Figure 5: (Color online) (a) Final spectrum of a linearly polarized input pulse (95% 
on the m — 1 axis and 5% along the m = 2 axis), as obtained by simulating the full 
Eq. (119[l . All other parameters for the pulse and the fiber are the same as in Fig. 
[4jc). Black solid and red dashed lines are respectively the polarization components of 
the input pulse along m — 1 and m = 2 respectively. Blue solid line is output pulse 
along m = 1, while green dashed line is along m = 2. (b) Time domain picture of 
the propagation shown in (a). The formation of several vector solitons of different 
amplitudes is observed, each of them subjected to RSFS according to their individual 
intensities. 
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